Conditional Modelling of Multivariate Extreme Value Data using R

Introduction
- Preliminaries
- Data
Exploratory multivariate modelling
- Exploratory plots
- Exploring pairwise extremal dependence
Conditional multivariate extreme value modelling
texmex 패키지를 활용한 conditional multivariate extreme value modeling
Constrained parameter space
Fitted model parameters
Prediction under the fitted model
- Building samples from multiple conditional models
- Joint exceedance curve estimation

Introduction

이 문서에서는 R 패키지 texmex를 가지고 다변량 자료의 극단값 분석을 어뗗게 하는지 보여준다. 이 분석 방법은 다음의 두 단계를 거친다.

Generalized Pareto distribution (GPD) modelling of the marginal variables followed by conditional multivariate extreme value modelling. 이 방법은 texmex1d 문서에 좀 더 자세히 설명되어 있다.

Preliminaries

texmex 패키지를 설치하면 library 커맨드를 이용해 패키를 불러온다.

library(texmex) 
#palette(c("black","purple","cyan","orange")) 
set.seed(20130618)
library(gridExtra)

Data

데이터는 Heffernan and Tawn (2004) 에 쓰인 자료를 사용한다. 이 자료는 다섯 가지 air pollutant들의 hourly means들의 daily maxima의 extremal behaviour를 묘사한다. 그 중 11월부터 2월까지의 자료를 모아놓은 winter 자료에 집중한다.

head(winter)

##   O3 NO2  NO SO2 PM10
## 1 27  50 112  13   34
## 2 27  51 126  13   29
## 3 15  43  90  21   33
## 4  9  71 470  44  101
## 5 20  51 167  48   30
## 6  8  50 211  16   44

summary(winter,digits=2)

##        O3          NO2            NO           SO2           PM10    
##  Min.   : 1   Min.   : 19   Min.   : 10   Min.   :  1   Min.   :  7  
##  1st Qu.:10   1st Qu.: 37   1st Qu.: 64   1st Qu.:  8   1st Qu.: 29  
##  Median :22   Median : 43   Median :112   Median : 15   Median : 40  
##  Mean   :20   Mean   : 44   Mean   :135   Mean   : 21   Mean   : 48  
##  3rd Qu.:29   3rd Qu.: 51   3rd Qu.:166   3rd Qu.: 26   3rd Qu.: 60  
##  Max.   :44   Max.   :104   Max.   :568   Max.   :200   Max.   :177

Exploratory multivariate modelling

당연히 일변량 극단값 분석보다 다변량 극단값 분석이 더 복잡하다. 일단 가장 먼저 떠오르는 이슈는 관찰값의 multivariate extreme을 어떻게 정의할 것이냐이다. 만약 관찰값이 모든 성분들에 대해 동시에 extreme을 갖는 경우만을 다변량 극단값으로 정의하면, 그러한 조건을 만족하는 자료는 충분치 않을 것이기 때문에 모델링을 하고 분석하기에 적합하지 않을 것이다. 한편 자료의 몸통(평균 근처) 부분을 가지고 dependency를 논하기에는 극단값에서의 dependence와 다를 것이기 때문에 항상 경향이 같다고 말할 수 없다.

Exploratory plots

먼저 우리는 자료의 변수들끼리의 pairwise dependence plot을 그려보기로 한다.

GGally::ggpairs(winter)

## Registered S3 method overwritten by 'GGally':
##   method from   
##   +.gg   ggplot2

이제 우리는 가로축에 NO를 넣고 나머지 변수들과의 산점도를 그려본다.

p1 <- ggplot(winter,aes(NO,O3)) + geom_point(colour="darkblue",alpha=0.5) 
p2 <- ggplot(winter,aes(NO,NO2)) + geom_point(colour="darkblue",alpha=0.5) 
p3 <- ggplot(winter,aes(NO,SO2)) + geom_point(colour="darkblue",alpha=0.5) 
p4 <- ggplot(winter,aes(NO,PM10)) + geom_point(colour="darkblue",alpha=0.5)

grid.arrange(p1,p2,p3,p4,ncol=2)

이 산점도로부터 우리는 NO과 다른 오염물질들간의 상관관계가 매우 다름을 알 수 있다. 오존의 경우는 NO가 클 때 음의 상관관계를 갖지만, \(\text{NO}_{2}\)와 \(\text{PM}_{10}\)의 경우는 양의 상관관계를 갖고 있음을 알 수 있다.

Exploring pairwise extremal dependence

다음은 Coles, Heffernan, and Tawn (1999) 에서 묘사된 summary statistics \(\chi\)와 \(\bar{\chi}\)를 이용해 pairwise extremal dependence를 체크해 보기로 한다.

chiO3 <- chi(winter[, c("O3", "NO")]) 
ggplot(chiO3, main=c("Chi"="Chi: O3 and NO", "ChiBar"="Chi-bar: O3 and NO"))

결과를 보면 두 번째 plot이 greyed out 되어 있는 것을 볼 수 있는데 이것은 quantile이 1로 갈 때 limiting value of \(\chi\)의 신뢰구간이 1을 초과한다는 것이다. 이는 asymptotic independence를 의미한다. This is evidence of asymptotic independence, in which case the plot of \(\chi\) is not relevant – since this shows the level of dependence only within the asymptotic dependence class.

chiNO2 <- chi(winter[, c("NO2", "NO")]) 
ggplot(chiNO2, main=c("Chi"="Chi: NO2 and NO", "ChiBar"="Chi-bar: NO2 and NO"))

이 그림은 다음과 같은 것들을 설명해준다.

Look at limiting value of \(\bar{\chi}(u)\) plot as the quantile \(u\) tends to 1: 1로 간다는 것은 asymptotic dependence의 증거다.
If limit in a. is equal to 1: asymptotic dependence 클래스 내에서의 dependence의 strength가 어느 정도나 되는지 살펴보기 위해 \(\chi(u)\)를 살펴본다. Quantile \(u\rightarrow 1\)일 때 이 함수의 극한값이 dependence의 strength를 의미하며 1에 가까울 수록 stronger dependence를 의미한다.
If limit in a. is less than 1: asymptotic dependence 클래스 내에서의 dependence의 strength가 어느 정도나 되는지 살펴보기 위해 \(\bar{\chi}(u)\)를 살펴본다. Although at asymptotic levels, the largest values of the variables tend not to occur in the same observation, at moderately extreme levels, dependence may still be relatively strong. Quantile \(u\rightarrow 1\)일 때 이 함수의 극한값은 dependence의 strength를 말해주는데 1에 가까운 양수일수록 positive dependence를, -1에 가까운 음수일수록 negative dependence를 의미한다. 0 근처 값은 near independence임을 의미한다.

위의 내용들을 앞선 그림에 적용해보면,

\(\text{O}_{3}\)와 NO의 경우 \(\bar{\chi}\)를 살펴보면 이들이 asymptotically independent이며 약한 negative dependence를 갖고 있음을 알 수 있다. 여기서는 \(\chi\) plot을 살펴보지 않아도 되며 texmex 패키지에서도 회색으로 처리한다.
\(\text{NO}_{2}\)와 NO의 경우 \(\bar{\chi}\) plot이 오른쪽으로 상향하며 1을 possible limit로 포함하고 있으므로 asymptotic dependence의 증거가 된다. The \(\chi\) plot indicates moderate positive dependence within this class.

Pairwise extremal dependence를 체크하는 다른 방법으로는 multivariate conditional Spearman’s correlation coefficient across a sliding window of values of the variables 가 있다.

mcsO3 <- MCS(winter[, c("O3", "NO")]) 
mcsNO2 <- MCS(winter[, c("NO2", "NO")])
g1 <- ggplot(mcsO3, main="MCS: O3 and NO") 
g2 <- ggplot(mcsNO2, main="MCS: NO2 and NO")
gridExtra::grid.arrange(g1,g2,ncol=1)

신뢰구간은 bootMCS를 이용해 추가할 수 있다.

bootmcsO3 <- bootMCS(winter[, c("O3", "NO")],trace=1000) 
bootmcsNO2 <- bootMCS(winter[, c("NO2", "NO")],trace=1000) 
g1 <- ggplot(bootmcsO3, main="MCS: O3 and NO")
g2 <- ggplot(bootmcsNO2, main="MCS: NO2 and NO")
gridExtra::grid.arrange(g1,g2,ncol=1)

multivariate conditional Spearman’s \(\rho\)의 plot이 같은 \(y\)축을 갖고 있지 않다는 것에 유의하며 해석하면, \(\chi\), \(\bar{\chi}\)와 같은 결과를 얻을 수 있다.

이 exploratory summaries들은 다음과 같은 사실을 얘기해준다. Conditional multivariate extreme value modeling을 할 때에 우리는 \(\text{O}_{3}\)와 NO 사이에는 negative association을 예상할 수 있고, \(\text{NO}_{2}\)와 \(\text{NO}\)를 분석할 때에는 stronger positive association을 예상할 수 있는 것이다.

Conditional multivariate extreme value modelling

Heffernan and Tawn (2004) 의 방법은

Fitting Generalised Pareto distribution (GPD) models to the marginal variables
Estimating the dependence structure

의 두 단계로 이루어진다.

Marginal transformation

The structure of the regression type dependence model is defined not on the original data scale, but after marginal transformation to standardised margins.

이 marginal transformation에서도 여러 가지 방법이 있다.

Heffernan and Tawn (2004) 는 Gumbel margin을 사용
그러나 Laplace margin을 사용하면 모수를 줄일 수 있다.

\(\mathbf{X}=(X_1, \ldots , X_d)\)를 \(d\)차원 확률변수이며 arbitrary marginal distribution을 갖는다고 하자.

\(\hat{F}_{i}\)를 \(i\) \((i=1,\ldots, d)\)번째 marginal distribution의 estimate라고 하고
\(G\)를 distribution function of the standardised marginal distribution 이라고 하자 (이것을 굼벨 또는 Laplace로 놓는 듯 하다).

Original vector variable \(\mathbf{X}\)가 \(\mathbf{Y}=(Y_1,\ldots, Y_d)\)로 바뀐다고 할 때 다음의 probability integral transform을 쓰면 된다.

\[\begin{equation} Y_{i} = (G^{-1}(\hat{F}_{i}(X_i)), \qquad{i=1, \ldots, d.} \tag{1} \end{equation}\]

실전에서는 \(\hat{F}_{i}\)를

자료의 empricial distribution function을 쓰거나
semi-parametric model using the empirical distributions below a threshold and the fitted GPD models for the tails of the distributions above the threshold

를 쓴다.

Regression model structure

\(Y_{i}, i\in \{1,\ldots, d\}\)를 the variable on which we are to condition 이라고 두자. 그러면 \(\mathbf{Y}_{-i}\)를 벡터 \(\mathbf{Y}\)에서 \(i\)번째 성분을 제외한 나머지 벡터라고 둘 수 있다. Heffernan and Tawn (2004) 의 방법은 \(Y_i\)가 어떤 high threshold \(t\)보다 높다고 conditioning을 한 후에, 즉 \(Y_{i}>t\)임을 관찰하였다고 한 후 남아있는 \(\mathbf{Y}_{-i}\)의 dependence를 모델링하는 것이다.

이런 conditional dependence structure의 회귀 모형의 형태는 식 (1) 에서 \(G\)를 어떻게 선택하느냐에 따라 달라진다.

Laplace margins: 이 경우는 \(G\)를 Laplace distribution function, \(\mathbf{Y}\)를 marginally Laplace distributed 되었다고 두는 것이다. \(Y_i\)가 high threshold \(t\)를 초과하였다고 conditioning을 하였을 때 남아있는 variable들 \(\mathbf{Y}_{-i}\)에 대한 Heffernan and Tawn (2004) 의 모형은 다음과 같은 형태를 갖는다.

\[\begin{equation} \mathbf{Y}_{-i} = \boldsymbol{\alpha}_{|i}Y_i + (Y_{i})^{\boldsymbol{\beta}_{|i}}\mathbf{Z}_{|i}, \tag{2} \end{equation}\]

이 때

\(\mathbf{Z}_{|i}\): residual vector
\(\boldsymbol{\alpha}_{|i}\): \((d-1)\) 차원 parameter vector, \(\boldsymbol{\alpha}_{|i} \in [-1,1]^{d-1}\). 참고로 \(0<\alpha_{j|i}\leq 1\)이면 \(Y_j\)와 large values of \(Y_i\) 사이에는 positive association을, \(-1\leq\alpha_{j|i}< 0\)일 때에는 negative association을 갖는다.
\(\boldsymbol{\beta}_{|i}\): \((d-1)\) 차원 parameter vector, \(\boldsymbol{\beta}_{|i} \in (-\infty, 1)^{d-1}\)

Gumbel margins: 이 경우는 \(G\)를 Gumbel distribution function, \(\mathbf{Y}\)를 marginally Gumbel distributed 되었다고 두는 것이다. \(Y_i\)가 high threshold \(t\)를 초과하였다고 conditioning을 하였을 때 남아있는 variable들 \(\mathbf{Y}_{-i}\)에 대한 Heffernan and Tawn (2004) 의 모형은 다음과 같은 형태를 갖는다.

\[\begin{equation} \mathbf{Y}_{-i} = \boldsymbol{\alpha}_{|i}Y_i + I_{\boldsymbol{\alpha}_{|i}=0, \boldsymbol{\beta}_{|i}<0}(\mathbf{c}_{|i} - d_{|i} \log Y_{i}) (Y_{i})^{\boldsymbol{\beta}_{|i}}\mathbf{Z}_{|i}, \tag{3} \end{equation}\]

이 때

\(\mathbf{Z}_{|i}\): residual vector
\(\boldsymbol{\alpha}_{|i}\): \((d-1)\) 차원 parameter vector, \(\boldsymbol{\alpha}_{|i} \in [0,1]^{d-1}\). Laplace 케이스와 비교해봤을 때 negative dependence를 표현하여려면 \(\alpha_{j|i}=0\)이 되어야 함을 알 수 있다.
\(\boldsymbol{\beta}_{|i}\): \((d-1)\) 차원 parameter vector, \(\boldsymbol{\beta}_{|i} \in (-\infty, 1)^{d-1}\)
\(\mathbf{c}_{|i}\): \((d-1)\) 차원 parameter vector, \(\mathbf{c}_{|i}\in (-\infty, \infty)^{d-1}\)
\(\mathbf{d}_{|i}\): \((d-1)\) 차원 parameter vector, \(\mathbf{d}_{|i}\in (0,1)^{d-1}\).

\(Y_{j}\)와 large \(Y_{i}\) 사이의 positive association을 나타내려면 \(\alpha_{j|i}>0\), \(\beta_{j|i}<)\)이어야 한다. Negative dependence는 \(\alpha_{j|i}=0\)이어야 하고 추가적으로 \(c_{j|i}, d_{j|i}\)에 대한 모델링을 해야 한다.

두 케이스를 비교해 봤을 때 Gumbel 대신 Laplace margin을 사용했을 때 모형이 훨씬 더 단순해짐을 알 수 있고 추론 또한 좀 더 직관적으로 (특히 weak dependence) 할 수 있음을 알 수 있다.

참고할 점은 두 경우 모두 model residuals \(\mathbf{Z}_{|i}\)에 대한 특별한 분포가정을 하지 않는다는점이다. 따라서 Heffernan and Tawn (2004) conditional dependence modeling을 semi-parametric이라고 할 수 있는 것이다.

즉 Heffernan and Tawn (2004) 의 방법론 하에서 conditioning variable \(Y_i\)와 remaining variables \(\mathbf{Y}_{-i}\) 사이의 에dependence를 완벽히 묘사하기 위해서는

parametric regression type model (식 (2) 또는 (3))
model residuals \(\mathbf{Z}_{|i}\)의 distribution modeling (Heffernan and Tawn (2004) 에서는 observed model residuals의 empirical distribution으로 처리했으나 후속 연구들에서는 다른 모델링을 적용하기도 한다)

These model residuals are calculated by using transformed data \(\mathbf{Y}\) and estimates of model parameters \(\hat{\boldsymbol{\alpha}}, \hat{\boldsymbol{\beta}}\) (and possibly also \(\hat{\mathbf{c}}\) and \(\hat{\mathbf{d}}\)) in (2) or (3) .

Constraints on parameter space

Heffernan and Tawn (2004) 의 후속연구인 Keef, Papastathopoulos, and Tawn (2013) 에서는 적합한 모형의 validity 이슈에 대해 다루었다. Keef, Papastathopoulos, and Tawn (2013) 논문의 내용은 fitted model이 유효하기 위해서는 위에서 묘사한 Heffernan and Tawn (2004) 의 constraint보다 좀 더 tight한 constraint가 필요하다는 것이다.

Keef, Papastathopoulos, and Tawn (2013) 에서 제안한 constraint는 fitted dependence model with the strength of extremal dependence exhibited by the data를 enforce한다. 이렇게 만들어진 constraint는 dependence parameter space의 boudary를 곡선으로 만들어준다. 다행히 texmex 패키지에서는 Laplace margin 하에서 이러한 constraint estimation을 제공한다. Diagnostic plot 또한 제공한다.

`texmex` 패키지를 활용한 conditional multivariate extreme value modeling

Model fitting

여기서는 Heffernan and Tawn (2004) 에도 쓰였던 winter 데이터셋을 쓰기로 한다. 아래 코드에서 mqu는 marginal GPD 모형을 적합할 때 threshold를 결정할 marginal quantile을 설정해 주는 것이다.

mex.O3 <- mex(winter, mqu=.7, penalty="none", which="O3")

## Warning in mexDependence(x = res1, which = which, dqu = dqu, margins = margins, : Assuming same quantile for dependence thesholding as was used
##      to fit corresponding marginal model...

mex.NO2 <- mex(winter, mqu=.7, penalty="none", which="NO2")

## Warning in mexDependence(x = res1, which = which, dqu = dqu, margins = margins, : Assuming same quantile for dependence thesholding as was used
##      to fit corresponding marginal model...

mex.NO <- mex(winter, mqu=.7, penalty="none", which="NO")

## Warning in mexDependence(x = res1, which = which, dqu = dqu, margins = margins, : Assuming same quantile for dependence thesholding as was used
##      to fit corresponding marginal model...

mex.SO2 <- mex(winter, mqu=.7, penalty="none", which="SO2")

## Warning in mexDependence(x = res1, which = which, dqu = dqu, margins = margins, : Assuming same quantile for dependence thesholding as was used
##      to fit corresponding marginal model...

mex.PM10 <- mex(winter, mqu=.7, penalty="none", which="PM10")

## Warning in mexDependence(x = res1, which = which, dqu = dqu, margins = margins, : Assuming same quantile for dependence thesholding as was used
##      to fit corresponding marginal model...

위에서 쓴 mex 함수는

migpd: marignal modeling
mexDependence: conditional dependence modeling

함수 두 개를 합쳐놓은 것이다. 즉 mex함수를 풀어서 쓰면 다음과 같이 되는 것이다. (conditioning variable이 O3일 경우)

marg <- migpd(winter, mqu=0.7, penalty="none")
mex.O3 <- mexDependence(marg, which = "O3")

## Warning in mexDependence(marg, which = "O3"): Assuming same quantile for dependence thesholding as was used
##      to fit corresponding marginal model...

사실 marginal modeling을 할 때의 threshold와 conditional dependence modeling을 할 때의 threshold를 결정할 quantile을 다르게 설정할 수 있는데, mexDependence에서 따로 입력을 하지 않으면 mqu를 그대로 쓰게 된다. 만약 다르게 하고 싶으면 다음과 같이 dqu를 따로 설정하면 된다 (아마 mqu<dqu여야 할 것이다).

mexDependence(marg, which = "O3", dqu=0.8)

## mexDependence(x = marg, which = "O3", dqu = 0.8)
## 
## 
## Marginal models:
## 
## Dependence model:
## 
## Conditioning on O3 variable.
## Thresholding quantiles for transformed data: dqu = 0.8
## Using laplace margins for dependence estimation.
## Constrained estimation of dependence parameters using v = 10 .
## Log-likelihood = -183.3875 -183.6389 -159.521 -168.4866 
## 
## Dependence structure parameter estimates:
##        NO2      NO     SO2     PM10
## a -0.07103 -0.1349 -0.2638 -0.12210
## b -0.46180 -0.3338 -0.1054 -0.08856

Marginal model diagnostics

이젠 적합된 marginal model의 diagnostics를 하는 방법에 대해 체크해보자. texmex 패키지에서는 mrlPlot, gpdRangeFit 함수를 제공한다.

g <- ggplot(marg)

do.call("grid.arrange", c(g[[1]], list(ncol=2, nrow=2))) # ... etc

do.call("grid.arrange", c(g[[2]], list(ncol=2, nrow=2))) # ... etc

ggplot(gpdRangeFit(winter$O3)) # ... etc

## Fitted values of xi < -0.5
## Fitted values of xi < -0.5
## Fitted values of xi < -0.5
## Fitted values of xi < -0.5
## Fitted values of xi < -0.5

## [[1]]

## 
## [[2]]

ggplot(mrl(winter$O3)) # ... etc

Dependence model diagnostics

O3 변수를 conditioning한 상황에서의 모델 결과를 불러보자.

ggplot(mex.O3)

## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'

이 그림은 위에서부터 차례로 (top) centered and scaled values of the dependence model residual across the range of extreme conditioning variable, (middle) absolute values of these, (bottom) and the original untransformed data with contours showing quantiles of the fitted conditional model을 보여준다.

만약 모형이 자료를 잘 적합한다명 첫 줄과 중간 줄의 그림에서는 어떠한 패턴도 나타나지 않아야 한다. 마지막 줄의 그림에서는 raw data distribution (세로줄 왼쪽) 과 fitted quantile의 shape가 잘 일치해야 한다.

다음은 NO 변수를 conditioning 했을 때의 결과이다.

ggplot(mex.NO)

## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'

대부분의 결과는 threshold 선택이 옳았음을 보여주나, 맨 윗쪽 SO2와 NO 사이의 그림을 보면 NO값이 커짐에 따라 centered and scaled된 residual이 감소하는 트렌드를 갑고 있을을 확인할 수 있다.

이번엔 threshold를 변화시켰을 때 dependence structure parameter estimate는 어떻게 변하는지에 대해 plotting 해 볼 수 있다. 적당히 높은 threshold에서 우리는 parameter estimate들이 constant이길 기대할 것이다. To gain some feeling for the variability in the parameters, we perform 10 (by default) bootstrap samples. trace=11로 놓은 것은 출력을 위해 임의로 바꾼 것이다. (The default is to report every ten replicates).

mrf <- mexRangeFit(marg, "NO", trace=11)

## Fitted values of xi < -0.5
## Fitted values of xi < -0.5
## Fitted values of xi < -0.5
## Fitted values of xi < -0.5
## Fitted values of xi < -0.5
## Fitted values of xi < -0.5
## Fitted values of xi < -0.5
## Fitted values of xi < -0.5
## Fitted values of xi < -0.5
## Fitted values of xi < -0.5
## Fitted values of xi < -0.5
## Fitted values of xi < -0.5
## Fitted values of xi < -0.5
## Fitted values of xi < -0.5
## Fitted values of xi < -0.5

ggplot(mrf)

이 그림에서는 NO2|NO 케이스 (두번째 줄) 에 문제가 있음을 확인할 수 있다. We try using different starting values. There is an option for using a previously fitted dependence model as a starting point, see the documentation for mexDependence.

start <- coef(mex.NO$dependence)[1:2,] # alternative starting value

mrf <- mexRangeFit(marg, "NO", trace=11,start=c(0.1,0.1))

## Fitted values of xi < -0.5
## Fitted values of xi < -0.5
## Fitted values of xi < -0.5
## Fitted values of xi < -0.5
## Fitted values of xi < -0.5
## Fitted values of xi < -0.5
## Fitted values of xi < -0.5
## Fitted values of xi < -0.5
## Fitted values of xi < -0.5
## Fitted values of xi < -0.5
## Fitted values of xi < -0.5
## Fitted values of xi < -0.5
## Fitted values of xi < -0.5
## Fitted values of xi < -0.5
## Fitted values of xi < -0.5
## Fitted values of xi < -0.5
## Fitted values of xi < -0.5
## Fitted values of xi < -0.5

ggplot(mrf)

새로 얻은 결과에서는 특별한 이상점을 찾을 수 없으며, 새로 starting value를 설정함으로써 문제를 해결할 수 있음을 보여준다. 즉 이 parameter estimate의 stability plot은 70퍼센트 percentile이 적절한 threshold임을 보여준다.

Constrained parameter space

우리가 적합한 모형을 예측에 사용하기 전에, 우리는 먼저 우리가 적합한 parameter estimate들이 추정에 사용된 objective function의 true maximum과 일치하는지 살펴봐야 한다. 이것은 optimization의 결과가 starting value의 영향르 받기 때문에 나타나는 문제이다. texmex에서는 추정한 parameter estimate들이 실제 maximum likelihood estimate에 수렴하는지 시각화로 보여줄 수 있다. (물론 저차원 한정, conditional spatial extreme case에서는?)

여기서는 NO가 주어졌을 때 NO2를 모델링 한 것에 대한 결과만 출력하기로 한다. Depedence model parameter의 추정에 관한 profile-likelihood surface를 보여주기 위해 mexDependence 함수의 PlotLikDo=TRUE 옵션을 활용할 수 있다.

par(mfrow=c(3,4), mar=par("mar")/2)
marg.NO2.NO <- migpd(winter[,c("NO2","NO")],mqu=0.7)
mex.NO2.NO <- mexDependence(marg.NO2.NO, which = "NO", dqu=0.7, PlotLikDo=TRUE)

이 그림의 결과는 point estimate가 permissible parameter space의 모서리에 위치함을 보여준다.

par(mfrow=c(1,1))
mex.NO2.NO <- mexDependence(marg.NO2.NO, which = "NO", dqu=0.7, PlotLikDo=TRUE, PlotLikRange=list(a=c(0.6,0.8),b=c(0.1,0.3) ))

구역을 제한해서 다시 그렸을 때에도 point estimate가 objective function의 maximum이 됨을 확인해 볼 수 있다. 물론 다양한 starting point를 활용하면서 최종 수렴값이 초깃값의 영향을 받는지 체크해 볼 수도 있다. (mexDependence의 documentation 참고)

mexDependence(marg,which="O3",dqu=0.7,PlotLikDo=TRUE)

## mexDependence(x = marg, which = "O3", dqu = 0.7, PlotLikDo = TRUE)
## 
## 
## Marginal models:
## 
## Dependence model:
## 
## Conditioning on O3 variable.
## Thresholding quantiles for transformed data: dqu = 0.7
## Using laplace margins for dependence estimation.
## Constrained estimation of dependence parameters using v = 10 .
## Log-likelihood = -257.702 -256.6681 -231.7684 -234.0916 
## 
## Dependence structure parameter estimates:
##       NO2       NO     SO2     PM10
## a 0.01301 -0.07278 -0.1683 -0.04719
## b 0.02020  0.03038 -0.1418  0.07142

Fitted model parameters

만족할 만한 모형 적합 결과를 얻었다면 추정된 모형의 모수값들을 체크해 볼 시간이다. Dependence structure의 모수들은 a값이 1에 가까우면 positive extremal dependence, -1에 가까우면 negative extremal dependence를 나타내지만 그 이상으로 직관적으로 설명할 수는 없다.

mex.O3

## mexDependence(x = marg, which = "O3")
## 
## 
## Marginal models:
## 
## Dependence model:
## 
## Conditioning on O3 variable.
## Thresholding quantiles for transformed data: dqu = 0.7
## Using laplace margins for dependence estimation.
## Constrained estimation of dependence parameters using v = 10 .
## Log-likelihood = -257.702 -256.6681 -231.7684 -234.0916 
## 
## Dependence structure parameter estimates:
##       NO2       NO     SO2     PM10
## a 0.01301 -0.07278 -0.1683 -0.04719
## b 0.02020  0.03038 -0.1418  0.07142

출력된 결과의 a값으로 보아 네 가지 변수들 중 SO2가 O3의 large value들에 대해 가장 negatively dependent이고 나머지 변수들은 O3의 large value들과 아주 약한 음/양 상관관계를 갖고 있음을 확인할 수 있다.

mex.NO

## mex(data = winter, which = "NO", mqu = 0.7, penalty = "none")
## 
## 
## Marginal models:
## 
## Dependence model:
## 
## Conditioning on NO variable.
## Thresholding quantiles for transformed data: dqu = 0.7011278
## Using laplace margins for dependence estimation.
## Constrained estimation of dependence parameters using v = 10 .
## Log-likelihood = -230.9793 -225.9886 -220.2708 -243.5154 
## 
## Dependence structure parameter estimates:
##        O3     NO2     SO2     PM10
## a -0.2558 0.08888  0.2852  0.70820
## b -0.2986 0.45710 -0.3102 -0.09794

다음은 NO를 conditioning 했을 때의 결과인데, NO2, SO2, PM10 모두 NO의 large value들과 positive extremal dependence를 갖고 있음을 확인할 수 있다. 셋 중에는 PM10이 NO의 large value들과 가장 강한 positive extremal dependence를 갖고 있음을 확인할 수 있다. O3는 NO의 large value들과 weak negative dependence를 갖고 있음을 확인할 수 있다.

Prediction under the fitted model

모수들의 짝 \((a,b)\)로 묘사된 variable들의 pair들의 dependence는 항상 \(\mathbf{Z}_{|i}\)의 empirical distribution과 연관이 있다.

texmex의 predict 함수에서는 fitted conditional multivariate model에 importance sampling을 적용해 에측을 한다. 여기서는 NO의 90퍼센트 percentile (pqu=0.9) 위로 conditioning 했을 때 값들을 simulate 한다. (pqu는 conditioning variable의 percentile인 듯)

set.seed(20130619)
nsim <- 1000
pO3 <- predict(mex.O3, pqu=.9, nsim=nsim) 
pNO2 <- predict(mex.NO2, pqu=.9, nsim=nsim) 
pNO <- predict(mex.NO, pqu=.9, nsim=nsim) 
pSO2 <- predict(mex.SO2, pqu=.9, nsim=nsim) 
pPM10 <- predict(mex.PM10, pqu=.9, nsim=nsim)

얻어진 conditional distribution의 결과는 다음과 같다.

summary(pO3)

## predict.mex(object = mex.O3, pqu = 0.9, nsim = nsim)
## 
## Conditioned on O3 being above its 90th percentile.
## 
## 
## Conditional Mean and Quantiles:
## 
##      O3|O3>Q90 NO2|O3>Q90 NO|O3>Q90 SO2|O3>Q90 PM10|O3>Q90
## mean      36.6       42.2      89.1       12.2        37.5
## 5%        33.8       26.0      18.0        3.0        19.0
## 50%       36.1       43.0      78.5       10.0        32.0
## 95%       41.2       59.0     184.0       26.7        85.0
## 
## Conditional probability of threshold exceedance:
## 
##  P(O3>28|O3>Q90) P(NO2>49|O3>Q90) P(NO>149|O3>Q90) P(SO2>23|O3>Q90)
##                1             0.27            0.148            0.089
##  P(PM10>53|O3>Q90)
##              0.107

마지막의 conditional probability part를 보면 O3>28|O3>Q90, … PM10>53|O3>Q90 등으로 되어있는데 여기서 28, 53 등의 threshold는 GPD 모형을 적합할 때의 marginal threshold들 (즉 앞서 mqu=0.7로 정의되었으므로 0.7 quantile에 해당)이다. 하지만 summary function에서 mth를 임의로 정함으로써 이 threshold value들을 바꾸어줄 수 있다.

summary(pO3,mth=c(39,40,100,10,40))

## predict.mex(object = mex.O3, pqu = 0.9, nsim = nsim)
## 
## Conditioned on O3 being above its 90th percentile.
## 
## 
## Conditional Mean and Quantiles:
## 
##      O3|O3>Q90 NO2|O3>Q90 NO|O3>Q90 SO2|O3>Q90 PM10|O3>Q90
## mean      36.6       42.2      89.1       12.2        37.5
## 5%        33.8       26.0      18.0        3.0        19.0
## 50%       36.1       43.0      78.5       10.0        32.0
## 95%       41.2       59.0     184.0       26.7        85.0
## 
## Conditional probability of threshold exceedance:
## 
##  P(O3>39|O3>Q90) P(NO2>40|O3>Q90) P(NO>100|O3>Q90) P(SO2>10|O3>Q90)
##            0.177             0.57            0.414            0.465
##  P(PM10>40|O3>Q90)
##              0.287

다음의 그림은 importance samping으로 적합한 conditional model을 시각화해준다.

ggplot(pO3)

그림에는 원래 자료 (회색 점)들과 threshold 위해서 importance sampling으로 얻어진 자료 (prediction) (오렌지, 파랑 점들)을 보여준다. 파랑색 점은 예측 시점에서 conditioning variable이 unconditioning variable보다 큰 경우 (quantile 기준), 오렌지 색 점은 unconditioning variable의 값이 conditioning variable보다 큰 경우 (quantile 기준)이다.

각각의 그림에서 solid curve는 두 개의 변수들의 marginal distribution들의 equal quantile들을 이어서 그린 선이다. 두 변수들이 완벽하게 dependent하면 예측값들이 그 선 위에 놓여있을 것이다.

위에서 그린 그림들은 O3를 conditioning한 결과인데 모든 변수들의 상관관계가 약하거나 또는 음의 상관관계를 보이고 있음을 확인할 수 있다. 이것을 NO를 conditioning하였을 때의 결과와 비교해 보면 흥미롭다.

ggplot(pNO)

여기서는 겨울철 PM10과 NO가 strong extremal dependence를 갖고 있음을 확인할 수 있다. 그리고 예측 결과를 보면 PM10의 marginal quantile이 높은 값들 (오렌지)과 NO의 marginal quantile이 높은 값들 (파랑)이 적절히 잘 섞여있음을 알 수 있다.

predict 함수에서 쓰인 importance sample들은 또한 conditioning variable의 treshold 위쪽의 특정 구간에 놓일 확률을 예측할 때 쓸 수도 있다. 또는 다변량 변수들이나 functional을 계산해 낼 수도 있다.

Building samples from multiple conditional models

때때로 우리는 single component만 큰 경우가 아닌 whole of the joint distribution of the multivariate random variable로부터 sampling을 해야 할 경우가 있다. 이 sampling approach는 분포의 tail들의 특정한 영역에 위치한 failure set에 event들이 위치할 확률을 계산하는 데 사용될 수 있다.

여기서는 각각의 변수들을 conditioning한 상태로 얻어진 conditional model들의 collection으로 정의된 whole of the modelled joint distribution으로부터 large Monte Carlo sample을 어떻게 구성할 수 있는지를 보여준다.

이렇게 conditional model들을 함께 모으는 과정은 우선 각각의 변수들을 conditioning한 상태로 모형을 먼저 돌린 후에 취하게 된다.

여기서는 marginal threshold를 0.7로, dependence quantile thresholdfmf 0.7로 둔 후 Heffernan and Tawn (2004) 의 예를 이용할 것이다.

mAll <- mexAll(winter,mqu=0.7,dqu=rep(0.7,5))

이제 collection of fitted model들로부터 Monte Carlo sample을 generate할 것이다. Heffernan and Tawn (2004) 에서처럼 random vector의 \(i\)번째 성분에 대해 conditioning을 하는 모형을 사용하여 이 \(i\)번째 성분이 모든 구성 요소 중 가장 큰 (샘플링 스케일로 측정된) 표본공간의 해당 부분에서 값을 시뮬레이션한다는 것이다. 다음과 같이 하면 된다.

원래 자료로부터 Monte Carlo sample을 전체 데이터셋에서 uniformly with replacement로 required number of obs를 뽑아 생성한다.
단계 1에서 얻은 Monte Carlo sample을 fitted GPD 모형을 활용해 Laplace scale로 transform한다.
Laplace scale에서 각 transformed data point에서 어떤 성분이 가장 큰지를 구분해난다. (since we are working on the common Laplace scale, this step calculates which component represents the highest marginal quantile)
단계 3에서 identify한 Monte carlo sample 중 어떤 것이 corresponding conditional dependence model threshold 위에 있는지 구분해낸다. (예를 들면, \(i\)번째 성분이 가장 큰 성분인 모든 점들에서는 we find which of these have ith component above the dependence threshold used to fit the conditional model given the ith component is above a given threshold)
각각의 conditioning variable에 대해 차례로, fitted conditional distribution으로부터 large independent sample을 생성해낸다. (conditional upon being above the associated dependence model fitting threshold) This is carried out on the original scale of the data.
Original scale of the data 에서 각각의 변수들에 대해 차례로, replace those values in our Monte Carlo sample from step 1 which are both above their conditional model threshold and for which the conditioning variable is the largest component (in Step 4) by a value generated from the appropriate conditional model (in Step 5).

The following plots highlight the selection of points in step 4. Points fulfilling the requirements of step 4 are shown by blue dots:

시뮬레이션은 texmex 패키지의 mexMonteCarlo 함수를 작동시키면 된다. 여기서는 원래 자료 (below the dependence thresholds)와 the collection of conditional models above each of the dependence thresholds로부터 5000개의 point들을 generate한다.

mexMC <- mexMonteCarlo(5000, mAll)

## Warning in qexp(p, lower.tail = lower.tail, log.p = log.p): NaNs produced

## Warning in qexp(p, lower.tail = lower.tail, log.p = log.p): NaNs produced

## Warning in qexp(p, lower.tail = lower.tail, log.p = log.p): NaNs produced

## Warning in qexp(p, lower.tail = lower.tail, log.p = log.p): NaNs produced

## Warning in qexp(p, lower.tail = lower.tail, log.p = log.p): NaNs produced

각각의 margin에 대해 dataset that were replaced by points sampled parametrically from the corresponding conditional tail model로부터 original sample의 point들의 갯수는 다음과 같다.

mexMC$nR

##   O3  NO2   NO  SO2 PM10 
## 1336  703  588  743  505

이것은 상당히 많은 양의 샘플들이 O3가 다른 margin들보다 가장 extreme한 component인 point들로 대체되었음을 보여준다. (O3 has around twice as many points replaced as any of the other margins). 이것은 O3가 클 때 O3와 다른 변수들 사이에 매우 약한/ 또는 약간 음의 상관관계를 갖는 fitted model들에 대응된다.

mAll$O3$dependence

## Conditioning on O3 variable.
## Thresholding quantiles for transformed data: dqu = 0.7
## Using laplace margins for dependence estimation.
## Constrained estimation of dependence parameters using v = 10 .
## Log-likelihood = -257.7152 -256.7165 -231.7644 -234.1655 
## 
## Dependence structure parameter estimates:
##       NO2       NO     SO2     PM10
## a 0.01444 -0.07238 -0.1678 -0.04756
## b 0.02122  0.03100 -0.1404  0.07132

이것의 결론은 O3가 클 때, 다른 변수들은 크지 않다는 것이다. 또한 texmex 패키지에서는 large Monte Carlo sample을 그리고 그것을 original dataset과 비교할 수 있다.

pairs(mexMC$MCsample)

그러나 required joint distribution으로부터 large sample들을 생성하기 위해 이 방법을 사용하는 것은 분명한 한계점들이 있다.

가장 근본적인 문제는 collection of conditional model들이 consistent estimator를 주지 않거나 심지어 잘 정의되지도 않는 joint distribution을 줄 수도 있다는 것이다. Heffernan and Tawn (2004) 의 방법은 전적으로 emprical한 방법이고 condtional model들을 추정하는 데 사용된 자료들의 underlying joint 분포들의 유효성에 의존한다. 우리는 sample underlying data로부터 얻어진 이러한 모형들이 underlying joint structure를 잘 반영하길 바라고 따라서 근사적으로 consistent distribution을 만들어주길 바라지만 실제로 반드시 그렇게 해주리라는 보장은 없다. Liu and Tawn (2014) 가 이 문제를 해결할 방법을 제시하긴 했지만 아직 texmex에 구현되지는 않았다.
적합한 threshold를 고르는 것의 중요성 또한 서로 다른 estimated model들을 묶을 때 중요하다. Marginal 그리고 dependence threshold는 empirical model로부터 (e.g. below the GPD or conditional model threshold)으로부터 parametric 모형 (above the respective threshold)의 transition이 부드럽도록 선택되어야 한다. Resulting semi-parametric model의 성분들 사이의 required continuity가 부족할 경우 연관된 Monte Carlo sample들의 결과도 좋지 않을 수 있다.

Joint exceedance curve estimation

일변량 setting에서 return level curve는 변수의 margival distribution이 어떻게 extrapolate되는지 보여준다고 할 수 있다. Return period는 return level과 동급 또는 그 이상의 이벤트들이 발생하기까지의 expected waiting time이라고 해석할 수 있다.

2차원에서는 이러한 단순한 curve가 존재하지 않는다. 예를 들어 200년을 return period라고 잡았을 때 i.i.d. 데이터의 경우 이것은 \(1/(365\times 200)\)의 확률(exceedance probability)과 같다. Univariate setting에서 우리는 fitted distribution with this exceedance probability의 quantile을 계산해 낼 수 있고 이것은 return level이 된다. 2차원 또는 그 이상인 경우에는 어떤 작은 확률이 주어졌을 때 both/all variable들이 일어날 많은 수의 joint event들이 존재하게 된다. 따라서 이 경우에는 exceedance probability \(p\)가 주어졌을 때 set of event들의 곡선을 묘사하게 되는데 이것을 joint exceedance curve라고 한다.

\[ \{ (x_{1,p}, \ldots, x_{d,p}) : P(X_1 > x_{1,p}, \ldots, X_{d}>x_{d,p})= p\}. \]

texmex 패키지에서는 2차원의 경우에 대해 구현되어 있다. 이 곡선은 다양한 방법으로 추정될 수 있다.

from the original data: for relatively non-extreme curves only, within the range of observations seen in the data
from a single fitted conditional model: conditioning on one variable only being large: curves can therefore be estimated only in that part of the space in which this estimated model is defined
from a collection of conditional models: each fitted conditional model is used to estimate that part of the joint exceedance curve in which that model’s conditioning variable is largest

여기서는 winter 데이터셋의 (NO2, NO) 변수를 활용해보기로 한다.

Joint exceedance curve directly from original data

아래 그림에서 나타나는 curve들은 joint exceedance probability들이 각각 0.2, 0.1, 그리고 0.05에 해당하는 경우이다. Raw 데이터들만 사용하는 경우에는 자료에서 관찰된 값 이상의 level들에 외삽할 수 없다. 이 곡선들은 empricial하게 추정되며 data의 숫자가 적어지는 high level에서는 결과가 좋지 않다.

WinterNO.NO2 <- winter[,3:2] #conditional이 아닌 joint인 듯
j1 <- JointExceedanceCurve(WinterNO.NO2,0.2) 
j2 <- JointExceedanceCurve(WinterNO.NO2,0.1) 
j3 <- JointExceedanceCurve(WinterNO.NO2,0.05)
ggplot(WinterNO.NO2,aes(NO,NO2)) + geom_point(colour="dark blue",alpha=0.5) + geom_jointExcCurve(j1,colour="orange") + geom_jointExcCurve(j2,colour="orange") + geom_jointExcCurve(j3,colour="orange")

Joint exceedance curve from fitted conditional model: one conditioning variable only

이번엔 fitted conditional model을 활용해 외삽해보고 우리가 곡선을 추정하고자 하는 joint tail region으로부터 importance sample을 생성해보고자 한다. 여기서는 extreme NO가 주어졌을 때 NO2의 conditional behavior를 계산해보기로 한다.

Importance sample들을 생설할 때 쓸 threshold (predict 함수의 pqu)을 유저가 알아서 설정할 수 있다.

p1 <- predict(mex.NO2.NO,nsim=5000,pqu=0.999) 
g <- ggplot(p1,plot.=FALSE)
j4 <- JointExceedanceCurve(p1,0.0005,which=c("NO","NO2"))
j5 <- JointExceedanceCurve(p1,0.0002,which=c("NO","NO2"))
j6 <- JointExceedanceCurve(p1,0.0001,which=c("NO","NO2"))
pl <- g[[1]] + geom_jointExcCurve(j4,aes(NO,NO2),col="purple") + geom_jointExcCurve(j5,aes(NO,NO2),col="purple") + geom_jointExcCurve(j6,aes(NO,NO2),col="purple")
pl

즉 앞 section과 달리 외삽을 하기 때문에 더욱 extreme level에 대한 joint exceedance curve를 그릴 수 있다.

또한 아래 그림처럼 그림 하나에 여러 개의 threshold 선택으로부터 얻은 importance sample들을 합쳐서 그릴 수도 있다.

p2 <- predict(mex.NO2.NO,nsim=5000,pqu=0.99)
j7 <- JointExceedanceCurve(p2,0.0005,which=c("NO","NO2")) 
j8 <- JointExceedanceCurve(p2,0.0002,which=c("NO","NO2")) 
j9 <- JointExceedanceCurve(p2,0.0001,which=c("NO","NO2"))
pl + geom_jointExcCurve(j7,aes(NO,NO2),col="purple") + geom_jointExcCurve(j8,aes(NO,NO2),col="purple") + geom_jointExcCurve(j9,aes(NO,NO2),col="purple")

계산된 joint exceedance curve들은 면시적으로 반환할 수 있는데, JointExceedanceCurve 함수의 x argument를 이용하면 된다.

Curve <- JointExceedanceCurve(p2,0.0005,which=c("NO","NO2"),x=seq(43,96,by=3)) 
Curve

## 
##  Estimated curve with constant joint exceedance probability equal to 5e-04 
##    NO      NO2
## 1  43 97.31809
## 2  46 97.31809
## 3  49 97.31809
## 4  52 97.31809
## 5  55 97.31809
## 6  58 97.31809
## 7  61 97.31809
## 8  64 97.31809
## 9  67 97.31809
## 10 70 97.31809
## 11 73 97.31809
## 12 76 97.31809
## 13 79 97.31809
## 14 82 97.31809
## 15 85 97.31809
## 16 88 97.31809
## 17 91 97.31809
## 18 94 97.31809

Joint exceedance curve from family of conditional models

필요하다면 여러 개의 변수들을 차례로 conditioning해서 적합한 여러 개의 모델을 묶을 수도 있다. 이러한 모형들의 집합의 fitting은 Section ?? 에서 다루었다. 우리는 Section ?? 에서 얻은 fitted model, 즉 object mAll을 sampling에 이용한다.

p3 <- mexMonteCarlo(nSample=5000,mexList=mAll)

## Warning in qexp(p, lower.tail = lower.tail, log.p = log.p): NaNs produced

## Warning in qexp(p, lower.tail = lower.tail, log.p = log.p): NaNs produced

## Warning in qexp(p, lower.tail = lower.tail, log.p = log.p): NaNs produced

## Warning in qexp(p, lower.tail = lower.tail, log.p = log.p): NaNs produced

## Warning in qexp(p, lower.tail = lower.tail, log.p = log.p): NaNs produced

j10 <- JointExceedanceCurve(p3,0.05,which=c("NO","NO2"))
j11 <- JointExceedanceCurve(p3,0.02,which=c("NO","NO2"))
j12 <- JointExceedanceCurve(p3,0.01,which=c("NO","NO2"))
ggplot(as.data.frame(p3$MCsample[,c("NO","NO2")]),aes(NO,NO2)) + geom_point(col="light blue",alpha=0.5) + geom_jointExcCurve(j10,aes(NO,NO2),col="orange") + geom_jointExcCurve(j11,aes(NO,NO2),col="orange") + geom_jointExcCurve(j12,aes(NO,NO2),col="orange")

## Warning: Removed 24 rows containing missing values (geom_point).

Smaller exceedance probability (아마도 extreme)에 대해서는 이러한 곡선들을 추정할 때 사용할 샘플의 크기를 임의로 크게 만들어 정확성을 높힌다.

마찬가지로 JointExceedanceCurve라는 것을 통해 curve를 만드는 점들의 값들을 반환할 수 있다. 앞선 코드에서 j10이 exceedance probability 0.05에 해당하는 curve를 그려주었으므로 아래 코드에서는 그에 해당하는 specific point value들을 반환한다.

JointExceedanceCurve(p3,0.05,which=c("NO","NO2"),x=seq(10,360,by=50))

## 
##  Estimated curve with constant joint exceedance probability equal to 0.05 
##    NO      NO2
## 1  10 66.31018
## 2  60 66.13116
## 3 110 66.12879
## 4 160 65.19954
## 5 210 64.28625
## 6 260 61.91513
## 7 310 58.59107
## 8 360 37.00000

Specifying return periods in terms of units of time

texmex 패키지에서 return period의 unit들은 observation이다. 왜냐면 몇몇 응용 사례에서 observation은 시간을 나타낼 수 있고, 다른 경우에는 예를 들면 개별 환자를 나타내는 경우에는 시간 스케일과는 전혀 관련이 없기도 하다.

그러나 몇몇 사례에서는 return level을 particular temporal return period들 (또는 joint exceedance curves to have a given Annual Exceedance Probability (AEP) of e.g. 1 in 10 years). 여기서는 200 year joint exceedance curve를 그려보기로 한다. 참고로 winter 데이터는 120년 관찰값이다.

ReturnPeriodInYears <- 200
NobsPerYear <- 120
ExceedanceProb <- 1/ (ReturnPeriodInYears * NobsPerYear) 
ExceedanceProb

## [1] 4.166667e-05

j200 <- JointExceedanceCurve(p1,ExceedanceProb,which=c("NO","NO2"), x=seq(700,by=50,len=5))
j200

## 
##  Estimated curve with constant joint exceedance probability equal to 4.166667e-05 
##    NO       NO2
## 1 700 110.36860
## 2 750 108.89374
## 3 800 105.63464
## 4 850  98.39811
## 5 900  81.86603

Plotting

j200 <- JointExceedanceCurve(p1,ExceedanceProb, which=c("NO","NO2"))
ggplot(WinterNO.NO2,aes(NO,NO2))+ geom_point(colour="light blue",alpha=0.5) + geom_jointExcCurve(j200,aes(NO,NO2),col="purple") + labs(title="200 year joint exceedance curve")

Coles, Stuart, Janet Heffernan, and Jonathan Tawn. 1999. “Dependence Measures for Extreme Value Analyses.” Extremes 2 (4): 339–65. https://doi.org/10.1023/a:1009963131610.

Heffernan, Janet E., and Jonathan A. Tawn. 2004. “A Conditional Approach for Multivariate Extreme Values (with Discussion).” Journal of the Royal Statistical Society: Series B (Statistical Methodology) 66 (3): 497–546. https://doi.org/10.1111/j.1467-9868.2004.02050.x.

Keef, Caroline, Ioannis Papastathopoulos, and Jonathan A. Tawn. 2013. “Estimation of the Conditional Distribution of a Multivariate Variable Given That One of Its Components Is Large: Additional Constraints for the Heffernan and Tawn Model.” Journal of Multivariate Analysis 115 (March): 396–404. https://doi.org/10.1016/j.jmva.2012.10.012.

Liu, Y., and J. A. Tawn. 2014. “Self-Consistent Estimation of Conditional Multivariate Extreme Value Distributions.” Journal of Multivariate Analysis 127 (May): 19–35. https://doi.org/10.1016/j.jmva.2014.02.003.

Conditional Modelling of Multivariate Extreme Value Data using R

Seoncheol Park

1/27/2020

Introduction

Preliminaries

Data

Exploratory multivariate modelling

Exploratory plots

Exploring pairwise extremal dependence

Conditional multivariate extreme value modelling

Marginal transformation

Regression model structure

Constraints on parameter space

`texmex` 패키지를 활용한 conditional multivariate extreme value modeling

Model fitting

Marginal model diagnostics

Dependence model diagnostics

Constrained parameter space

Fitted model parameters

Prediction under the fitted model

Building samples from multiple conditional models

Joint exceedance curve estimation

Joint exceedance curve directly from original data

Joint exceedance curve from fitted conditional model: one conditioning variable only

Joint exceedance curve from family of conditional models

Specifying return periods in terms of units of time

Conditional Modelling of Multivariate Extreme Value Data using R

Seoncheol Park

1/27/2020

Introduction

Preliminaries

Data

Exploratory multivariate modelling

Exploratory plots

Exploring pairwise extremal dependence

Conditional multivariate extreme value modelling

Marginal transformation

Regression model structure

Constraints on parameter space

texmex 패키지를 활용한 conditional multivariate extreme value modeling

Model fitting

Marginal model diagnostics

Dependence model diagnostics

Constrained parameter space

Fitted model parameters

Prediction under the fitted model

Building samples from multiple conditional models

Joint exceedance curve estimation

Joint exceedance curve directly from original data

Joint exceedance curve from fitted conditional model: one conditioning variable only

Joint exceedance curve from family of conditional models

Specifying return periods in terms of units of time

`texmex` 패키지를 활용한 conditional multivariate extreme value modeling